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SUMMARY 


i 


The  numerical  solution  of  continuum  problems  in  unbounded 
regions  involves  two  essential  approximations:  first,  the 

continuum  must  be  approximated  by  a discrete  unit;  and  second, 
the  unbounded  domain  must  be  approximated  by  a finite  domain. 

The  first  problem  is  the  one  usually  studied  in  numerical  analy- 
sis. The  second  has  received  much  less  attention  and  is  the 
subject  of  the  present  paper. 

The  present  work  was  motivated  by  the  problem  of  the  numer- 
ical simulation  of  boundary  layer  flows  in  transition  and  turbu- 
lent regimes.  The  prototype  of  such  flows  is  the  flow  over  a 
semi-infinite  flat  plate  undergoing  transition  to  turbulence. 

The  geometry  of  this  three-dimensional  flow  is  infinite  in  three 
directions.  The  formulation  of  satisfactory  boundary  conditions 
is  simplest  in  the  cross-stream  directions.  On  both  theoretical 
and  experimental  grounds,  periodic  boundary  conditions  can  be 
justified  in  this  direction.  On  the  other  hand,  treatment  of 
the  downstream  direction  is  not  so  simple.  The  mapping  techniques 
of  this  paper  cannot  be  used  effectively  for  this  aspect  of  the 
transition  problem.  However,  the  techniques  developed  here  are 
appropriate  for  the  treatment  of  the  boundary  conditions  in  the 
direction  normal  to  the  boundary  layer.  j 

The  present  study  is  restricted  to  one  special  technique 
for  the  treatment  of  the  point  at  infinity:  coordinate  trans- 

formation of  the  infinite  domain  onto  a finite  region.  The 
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utility  of  mapping  methods  are  examined  for  six  model  problems. 


1 


two  of  which  are  critical  components  of  the  boundary-layer 


transition  study.  These  model  problems  are  (1)  the  solution 


of  the  one-dimensional  diffusion  equation  in  a semi-infinite 


region;  (2)  the  eigenvalues  of  the  quantum-mechanical  oscillator; 


(3)  the  eigenvalues  of  the  Orr- Sommer f eld  equations  for  the 


Blasius  boundary  layer;  (4)  the  calculation  of  the  Palkner-Skan 


boundary-layer  profiles;  (5)  solutions  of  the  wave  equation;  and 


(6)  shock-wave  solutions  to  Burger's  equation.  The  utility  of 


these  mappings  is  determined  by  comparing  the  numerical  solu- 


tions of  these  six  problems,  using  different  mappings,  with  the 


exact  solutions. 


It  is  concluded  that  mappings  are  an  effective  way  to  solve 


problems  in  infinite  domains  provided  that  the  solution  is  simple 


at  infinity.  If  the  solution  oscillates  at  infinity,  then  infinity 


must  be  an  essential  singularity  of  the  solution  and  mappings  fail. 


When  mapping  is  applicable,  the  proper  choice  of  mapping  should 


be  based  on  the  criterion  that  the  solution  to  the  problem  be 


smooth  in  the  mapped  coordinate. 


Mapping  techniques,  and  in  particular  algebraic  mappings. 


are  successful  in  the  simulation  of  boundary  layer  flows  because 


the  boundary  condition  in  the  direction  normal  to  the  boundary 


layer  is  that  the  flow  is  a simple  laminar  free  stream.  The 


algebraic  mapping  allows  substantial  economy  in  the  number  of 


mesh  points  used  in  obtaining  the  numerical  solution.  The  effi- 


ciency of  this  mapping  is  most  important  in  the  context  of  multi- 


dimensional numerical  hydrodynamics,  where  economy  in  resolution 


is  vital. 


J 
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1 . Introduction 


The  numerical  solution  of  continuum  problems  in  unbounded 
regions  involves  two  essential  approximations:  first,  the  con- 
tinuum must  be  approximated  by  a discrete  set;  and,  second,  the 
unbounded  domain  must  be  approximated  by  a finite  domain.  The 
first  problem  is  the  one  usually  studied  in  numerical  analysis. 

The  second  has  received  much  less  attention  and  is  the  subject 
of  the  present  paper.  We  restrict  the  present  study  to  one 
special  technique  for  the  treatment  of  the  point  of  infinity: 
coordinate  transformation  of  the  infinite  domain  into  a finite 
region.  One  of  the  principle  conclusions  of  the  paper  is  that, 
while  transformations  are  not  universally  valid,  there  is  a class 
of  problems  for  which  it  is  a very  useful  technique. 

The  present  work  was  motivated  by  the  problem  of  the 
numerical  simulation  of  boundary  layer  flows  in  transition  and 
turbulent  regimes.  The  prototype  of  such  flows  is  the  flow 
over  a semi- infinite  flat  plate  undergoing  transition  to  tur- 
bulence. The  geometry  of  this  three-dimensional  flow  (see 
Fig.  1)  is  infinite  in  three  directions.  The  formulation  of 
satisfactory  boundary  conditions  is  simplest  in  the  y direction. 

On  both  theoretical  and  experimental  grounds,  periodic  boundary 
conditions  can  be  justified  in  y . On  the  other  hand,  treatment 
of  the  downstream  x direction  is  not  so  simple.  The  mapping 
techniques  of  the  present  paper  cannot  be  used  effectively  for 
this  aspect  of  the  transition  problem.  Techniques  for  the  imposition 
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of  inflow  and  outflow  boundary  conditions  (which  are  appropriate 
in  the  x direction)  will  be  discussed  elsewhere.  However, 
the  techniques  developed  here  are  appropriate  for  the  treatment 
of  the  z direction  (normal  to  the  boundary  layer) . We  shall 
find  that  mapping  techniques  are  successful  in  z because  the 
boundary  condition  at  z = 08  is  that  the  flow  is  a 
simple  laminar  free  stream  (in  the  present  case,  uniform  flow) . 

The  idea  of  mapping  an  infinite  geometry  into  a finite 
one  is  not  original.  For  example.  Van  de  Vooren  and  Dijkstra 
[1]  successfully  applied  coordinate  transformations  to  the 
numerical  solution  of  laminar  flow  past  a flat  plate;  Davis  [2] 
applied  similar  techniques  to  laminar  flow  past  a parabola. 

We  will  examine  the  utility  of  mapping  methods  for  six 
model  problems,  two  of  which  are  critical  components  of  the 
boundary- layer  transition  study.  In  Sec.  2 we  study  the  solution 
of  the  one-dimensional  diffusion  equation  in  a semi-infinite 
region.  In  Sec.  3,  the  eigenvalues  of  the  quantum-mechanical 
harmonic  oscillator  are  found  using  mappings.  In  Sec.  4,  the 
eigenvalues  of  the  Orr-Sommerfeld  equation  for  the  Blasius  boundary 
layer  are  calculated,  while  in  Sec.  5 mapping  techniques  are 
applied  to  the  calculation  of  Falkner-Skan  boundary- layer  profiles. 
The  examples  of  Sects.  6 and  7 illustrate  the  limitations  of  map- 
ping techniques.  Finally,  we  summarize  some  heuristic  rules 
for  the  applicability  of  mappings. 

2.  One-Dimensional  Diffusion  Equation  in  a Semi-Infinite  Domain 

Consider  the  mixed  initial-boundary  value  problem 
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ut  uxx 

(2.1a) 

u(x,t)  = 

0 

ft 

1 A 

0 

(2.1b) 

u(0,t)  = 

sin  t t > 

0 

(2.1c) 

u(x,t) 

bounded  as 

x -*  °° 

(2. Id) 

One  particular  physical  realization  of  these  equations  is 
the  Rayleigh  shear  flow  in  the  neighborhood  of  an  oscillating 
flat  plate  [3].  As  t <*>  , the  exact  solution  to  (2.1)  is 

asymptotically 


[which  do  yield  conditional  stability  restrictions  of  the  usual 
2 

kind  At  = 0 (h  )]  lead  to  an  unclosed  set  of  equations.  For 
example,  the  second-order  semi-discrete  scheme 


9u (x, t)  _ u (x+h, t) -2u (x, t) +u (x-h , t) 


(2.4) 


involves  u(x+h,t)  for  every  x , so  that  a finite  number  of 
eauations  in  the  same  number  of  unknowns  is  never  obtained  in 
a finite  x-interval. 

The  most  obvious  way  to  avoid  the  latter  problem  is  to 
impose  a boundary  condition  at  an  artificial  boundary  x = L , 
like 


u (L , t)  = 0 . 


(2.5) 


If  L is  fixed,  the  solution  to  (2.4-5)  does  not  converge  as 
h -►  0+  to  the  solution  to  (2.1).  However  in  the  double  limit 
L -*•  + 00  , h -*•  0+  , convergence  is  achieved. 

Another  way  to  handle  this  problem  is  to  use  a non-uniform 
grid.  Such  a grid  is  obtained  by  first  mapping  the  semi-infinite 
region  0 < x < <*>  onto  the  finite  region  0 <_  2 < 1 and  then 


using  the  uniform  grid 


Zj  = j/J 


0 < j < J 


fS.6) 


The  boundary  condition  (2. Id)  becomes  simply 
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t 


■where  L is  a constant  scale  factor.  In  Fig.  2 we  plot  z 
versus  x for  the  exponential  map  (2.7)  with  various  values  of 
L . In  Fig.  3 we  give  similar  plots  for  the  algebraic  map  (2.8). 
The  points  on  the  curves  in  Figs.  2 and  3 indicate  the  values 
of  x with  Zj  = .04j  (J  = 25).  For  both  maps,  the  equivalent 
mesh  in  x is  nonuniform  with  the  most  rapid  variation 
occuring  with  x >>  L . 

The  exponential  map  (2.7)  gives  slightly  better  resolution 
near  x = 0 than  the  algebraic  map  (2.8)  while  the  algebraic  map 
gives  much  better  resolution  than  the  exponential  map  as  x -*■+«>  . 
In  fact, 


x+L 


< 1 - e 


-x/L 


for  all  x > 0 . Thus,  if  a uniform  grid  is  used  in  the  mapped 
variables  and  L is  fixed,  the  grid  point  z = Az  lies  at  a 
slightly  smaller  value  of  x with  the  exponential  map  than 


. — . 


rtf**  > 


with  the  algebraic  map.  Conversely,  the  grid  point  z = 1 - Az 
lies  at  much  larger  x with  the  algebraic  map  (x  ~ L/Az) 
than  with  the  exponential  map  (x  . L tn  1/Az) 

The  maps  (2.7)  and  (2.8)  are  especially  convenient  because 
they  yield  simple  expressions  for  derivatives.  With  the  ex- 
ponential map  (2.7),  derivatives  with  respect  to  x become 


3u(x,t)  _ 1-z  3u  (z, t) 
3x  L 3 z 


(2.9a) 


3 u(x, t)  1 fi  t 3 #i  \ 3u(z,t) 
— =5 T (1-z)  -gj  (1-z)  ^ 


3x 


3z 


(2.9b) 


With  the  algebraic  map  (8) , derivatives  with  respect  to  x 
become 


3u(x,t)  = (1-z)2  3u(z,t) 

3x  L 3z 


(2.10a) 


3 u (x,  t) 

a? 


(1-Z) 


_3_ 

3z 


(1-z) 


3u (z, t) 
3z 


(2.10b) 


With  these  transformations  and  centered  space  differencing  of 
all  z derivatives,  a closed  set  of  equations  is  obtained  for  the 

numerical  solution  of  problem  (2.1). 

Let  us  now  compare  some  numerical  results  for  the  solution 
of  (2.1)  using  the  three  methods  discussed  above:  restriction 

to  x <_  L , exponential  mapping  (2.7),  and  algebraic  mapping  (2.8). 
In  general  there  are  two  kinds  of  numerical  error:  truncation 


6 


error  and  mapping  error.  In  the  case  of  the  solution  using  the 

restricted  x-interval  0 <_  x < L , the  truncation  error  is  of 
2 

order  h and  the  mapping  error  due  to  neglect  of  the  interval 
L <_  x < “ and  imposition  of  u(L,t)  =0  is  of  order  exp(-L//J) 
[see  (2.2)].  It  turns  out  that,  in  order  to  obtain  an  absolute 

-3 

error  of  10  in  the  asymptotic  solution  (2.2)  for  t -*■  » at 
x = 1 , numerical  solutions  using  the  restricted  domain  method 
require  L > 10  and  h < 1/5  , or  a total  of  at  least  50  grid 
points  (see  Table  1) . 

The  algebraic  map  allows  accurate  results  to-  be  obtained  much 
more  efficiently.  For  example,  to  achieve  an  error  of  less  than 
10  3 as  t -*■  00  at  x = 1 requires  less  than  15  grid  points  using 
the  algebraic  map  (2.8)  with  L =.l  (see  Table  1). 

The  results  given  in  Table  1 show  that  the  algebraic  map  gives 
a much  better  representation  of  the  solution  for  large  x than  the 
results  obtained  by  restricting  the  domain  or  the  exponential  map. 
The  reason  for  this  behavior  is  simply  that  the  amplitude  of  the 
solution  (2.1)  is 

f,  (z)  = e— WI^U-zH 


in  terms  of  the  algebraic  map  (2.8)  while  it  is 


(z) 


(1-z) 


L//2 


' 


7 


1 1 


, 


in  terms  of  the  exponential  map  (2.7).  The  function  f^z)  and  all 


its  derivatives  vanish  at  z = 1 . On  the  other  hand,  f2'(l)  = “ 


(for  L = 1)  which  induces  a relatively  large  error  when  the  exponen- 


tial map  is  used.  Observe  from  Table  1 that  the  errors  obtained 


using  the  exponential  map  do  not  decrease  as  n as  h -*■  0+  when 


L < 2/2  . When  L = 2/1  , f2(z)  and  all  its  derivatives  are  finite 


»~4 


at  z = 1 ; the  error  is  less  than  10  as  t -*■  00  at  x = 1 with 
only  about  35  points.  In  general,  if  the  exact  solution  or  one  of 
its  z-derivatives  is  singular  at  z = 1 (or  x = °°)  , large 


numerical  errors  result. 


3.  Eigenvalues  of  the  quantum-mechanical  harmonic  oscillator 
The  eigenvalues  X of  the  Hermite  equation 


12 


u"  - jx  u = -Xu,  u bounded  as  |x|  -*•<*>  (3.1) 


are  X = n + ■=■  (n=0,l,2, . . . ) . The  corresponding  eigenfunctions 


12, 


are  the  Hermite  functions  exp(-±x  )Hen(x)  , where  HeQ(x)  = 1 , 


He^(x)  = x , He2(x)  = x -1  , and  so  on.  If  n is  even  the 


eigenfunctions  are  even  functions  of  x ; if  n is  odd,  the 
eigenfunctions  are  odd  functions.  We  shall  only  study  the  even 
eigenvalues  and  eigenfunctions. 

The  numerical  solution  of  (3.1)  requires  a method  for  handling 
the  boundary  conditions  at  ± °°  . Here  we  compare  two  methods 


applied  to  the  determination  of  the  even  modes.  In  the  first 

2 


method,  we  assume  that  the  function  u is  a function  of  x 
alone  and  require 


u (L)  = 0 


(3.2a) 


I 


8 


E]  I 

for  some  large  L . In  the  second  method,  we  make  the  algebraic  map 


z = 2 —7T — j “ 1 (3.2b) 

x^+LZ 


and  seek  the  solution  as  a bounded  function  of  z for  -1  < z < 1. 


For  both  methods  of  handling  the  boundary  conditions  at 
« , we  use  Chebyshev  series  [4]  to  represent  the  eigenfunction 
u(x)  . In  the  first  method,  u(x)  is  represented  as 


N 

u(x)  = l anT2n(x/L)  ' (3.3) 

n=0 


and  the  boundary  condition  (3.2a)  is  applied;  in  the  second 
method,  u(x)  is  represented  in  terms  of  the  mapped  variable 
(3.2b)  as 

N 

u (x)  = l anTn (z)  • (3.4) 

n=0 

1 


Here  TR(y)  is  the  Chebyshev  polynomial  of  degree  n defined 
by 


T (cos  0)  = cos  n0 
n 


The  details  of  the  application  of  Chebyshev  series  to  the 
numerical  solution  of  ordinary  and  partial  differential  equations 
are  given  in  [4-5] . 


i 
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In  Tables  2-3,  we  present  numerical  results  obtained  by  the 

methods  described  above  for  the  even  modes  having  exact  eigenvalues 
15  9 

2 * 1 * \ ' resPectivelY-  should  be  apparent  that  the  algebraic 
mapping  achieves  high  accuracy  much  more  efficiently  than  does 
simple  truncation  of  the  infinite  interval.  Notice  8that  for  a 
given  number  of  Chebyshev  polynomials  N there  is  an  optimal 
choice  of  scale  L that  gives  the  most  accurate  result.  Also 
observe  that  for  fixed  L there  is  rapid  (faster  them  algebraic) 
convergence  of  the  eigenvalues  as  N -*•  • in  both  methods.  This 
is  a general  property  of  Chebyshev  expansions.  However,  when 
truncation  to  | x | £ L is  used  the  eigenvalues  converge  to  the 
wrong  answer  unless  the  simultaneous  limit  L -*■  « is  also  taken. 

In  Pig.  4,  we  plot  the  eigenfunction  corresponding  to  the 
eigenvalue  -|  obtained  using  the  algebraic  mapping  (3.3)  for 
various  numbers  of  Chebyshev  polynomials  N . The  eigenfunctions 
are  all  normalized  by  u(0)  - -1  . Here  the  exact  eigenfunction 
of  (3.1)  with  E *»  ^ - (1-x2) exp (-jx2)  . Notice  the  very  rapid 

convergence  to  the  exact  eigenfunction  as  N increases. 

4.  Orr-Sommerfeld  Equation  for  Blasius  Flow 

The  Orr-Sommerfeld  equation  governs  two-dimensional  linear 
disturbances  to  incompressible  parallel  shear  flows.  We  assume 
that  the  (dimensionless)  undisturbed  flow  velocity  is  CJ(z)x 
where  x is  a unit  vector  in  the  x-direction  and  that  the 
z-component  of  the  perturbation  velocity  is  proportional  to  the 
real  part  of 


10 


OTKOSanw 


w(z)  e 


ia (x-ct) 


where  a , the  longitudinal  wavenumber,  is  assumed  real 
and  c (usually  complex)  is  the  phase  speed  of  the  disturbance 
propagating  in  x . If  Im(c)  > 0 , the  disturbance  grows 
in  time  and  the  flow  is  unstable.  It  may  be  shown  that  the 
linearized,  incompressible,  Navier-Stokes  equations  can  be 
reduced  to  the  Orr-Sommerfeld  equation 


2 2 2 
“ a2^  w = iaR[  (U(z)-c)^— ^ - a2^w  - U"  (z)wl 


(4.1) 


where  R is  the  Reynolds  number.  On  rigid  no-slip  walls,  the 
perturbation  velocity  must  satisfy 


w = w'  =0  . 


(4.2) 


The  problem  (4.1)  with  the  homogeneous  boundary  conditions  (4.2) 
is  an  eigenvalue  problem  for  the  phase  speed  c , assuming  a 
is  given. 

The  laminar  flow  over  a flat  plate  z = 0 , as  in  Fig.  1, 
satisfies  the  conditions  of  the  preceding  paragraph  except  for 
a slow  variation  in  x which  we  neglect;  U(z)  is  the  Blasius 
velocity  profile  which  is  determined  by  equations  summarized  in 
Sec.  5.  Linearized  disturbances  to  the  Blasius  flow  are  governed 
by  (4.1)  with  the  boundary  conditions  (4.2).  As  z -*■  °°  , the 
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disturbance  should  remain  bounded;  it  may  be  shown  that  this 
condition  becomes 


w(z)  ~ e-ctz  (z  -*■  «)  (4.3) 

if  R is  large  and  c is  not  close  to  1 . 

The  problem  is  to  solve  numerically  the  eigenvalue  problem 
(4.1)  - (4.3)  in  the  region  0 <_  z < ® . We  will  compare  several 

methods  for  handling  the  boundary  condition  at  « . First, 

the  region  may  be  truncated  to  the  region  0 £ z < L with  the 
artificial  boundary  conditions 

w (L)  = w'  (L)  = 0 (4.4) 

imposed  on  the  finite  lid  z = L . Second,  the  asymptotic 
behavior  (4.3)  may  be  used  to  infer  the  improved  boundary  condition 

w' (L)  + aw(L)  = 0 (4.5) 

on  z ■ L . The  other  methods  we  will  compare  use  the  exponential 
and  algebraic  mappings  (2.7-8)  of  the  semi-infinite  region 
0 <_  z < » into  a finite  domain. 

The  exponential  map 

Z = 1 - 2e_z/L  (4.6) 


— ----- 


types  of  boundary  conditions  will  be  applied  at  Z ■ 1 (z  = 08 ) 


Three  kinds  of 


boundary  conditions  will  be  applied  at  Z * 1 


We  have  computed  the  eigenvalues  of  the  Orr-Sommerfeld 


equation  for  these  various  cases  by  expansion  of  w in  series 
of  Chebyshev  polynomials  [4] , the  eigenvalues  being  found  by 
matrix  eigenvalue  methods  [6] . Comparisons  between  the  methods 
are  given  for  the  case  R = 580,  a = 0.179,  a case  previously 
studied  in  some  detail  [7,8]  . Results  are  given  in  Table  4 


only  for  the  single  unstable  eigenvalue  (Im  c > 0)  for 


this  choice  of  R , a . Here  N is  the  number  of  Chebyshev 

I 

polynomials  used  to  represent  the  modes.  High  resolution 

1 E 

numerical  calculations  yield  the  value 

c = 0.36412286  + i 0.00795972  (4.13) 


to  eight  decimal  places.  In  all  cases,  w(0)  = w' (0)  =0  at 
the  rigid  wall  z * Z » 0 . 

Comparison  of  cases  1-4  with  the  'exact'  solution  (4.13) 
shows  that  the  error  in  c incurred  by  truncation  to  z <_  L 
is  of  order  e”  ; also,  for  fixed  L , the  convergence 
with  increasing  N is  very  rapid,  consistent  with  the  expected 
infinite  order  rate  of  convergence  of  Chebyshev  series  [4] . Cases 
5,6  show  that  the  'improved'  boundary  condition  (4.5)  does  not 
help,  at  least  for  the  values  of  N used  in  the  present  cal- 
culations. In  fact,  use  of  (4.5)  does  improve  the  results  if  N 
is  large.  We  do  not  enter  into  these  comparisons  further  because 
they  are  not  central  to  the  present  paper. 

The  results  of  cases  7-10  for  the  exponential  map  are  very 
disappointing.  With  the  upper  boundary  condition  (4.7)  the 


f convergence  rate  is  roughly  like  1/N  ; the  boundary  conditions 


E 


Since  aL  = 0.179  for  cases  7-10,  it  follows  that  w has 

a cusp  singularity  at  Z = 1 . This  singularity  in  the 

derivative  of  w at  Z = 1 slows  the  convergence  of  the 

Chebyshev  series  representation  of  w . The  boundary  conditions 

(4.8)  are  inconsistent  with  the  asymptotic  behavior  (4.3)  since 

lim  w' (Z)  = ® ; nevertheless , convergence  to  the  correct  eigen- 

Z-*l- 

value  seems  to  be  obtained,  but  th^Tnconsistency  destroys  the 
rapid  convergence  properties  of  Chebyshev  series. 

On  the  other  hand,  the  algebraic  map  (4.9)  works  remarkably 
well.  Comparison  of  cases  11-16  with  the  'exact'  result  (4.13) 
shows  that  extremely  rapid  convergence  is  achieved  at  low  values 
of  N . The  reason  is  that  in  the  algebraically  mapped  co- 


ordinate (4.3)  becomes 


w ~ e 


-aL (1+Z) / (1-Z) 


(Z  - 1-) 


Thus,  w and  all  its  derivatives  with  respect  to  Z approach 
0 as  Z -*■  1—  . In  contrast  to  the  exponential  map,  the 
algebraic  map  gives  nicely  behaved  solutions  as  Z 1-  . 

5.  Application  to  the  Falkner-Skan  equation 

Consider  the  semi- inf inite  flat  plate  shown  in  Fig.  1 
in  which  the  uniform  inflow  at  x = 0 is  tilted  at  an  angle 
^ 8ir  with  respect  to  the  plate.  It  can  be  shown  that  the  x 
component  of  the  velocity  at  z = ® is 

U(x)  = U0(x/£)m  (5.1) 


' -*  ''  VSM*- ' * *•  X 


Here  UQ  is  the  free  stream  velocity  at  x = 0 , i is  a length 
scale  and 

m = 6/(2-0)  . (5.2) 

If  $ > 0 the  flow  is  accelerated  along  the  plate;  this 
case  can  also  be  interpreted  as  the  flow  over  a wedge  with  an 
included  angle  of  0tt  . If  8 < 0 the  flow  is  decelerated  and 
the  flow  is  just  the  same  as  occurs  on  the  underside  of  the 
plate  when  0 > 0 . Finally  if  0=0  the  flow  is  parallel 
to  the  plate  and  neither  accelerates  nor  decelerates. 

The  requirement  that  the  velocity  at  the  surface  of 
the  plate  be  zero  leads  to  the  formation  of  a boundary  layer  in 
the  immediate  vicinity  of  the  plate.  This  boundary  layer  flow 
problem  is  solved  by  introducing  the  similarity  variable 


n = z[|(m+l)o  (x)/vx]1/2  (5.3) 


and  a stream  function  ip(x,z)  related  to  the  x and  z com- 
ponents of  the  velocity,  u and  v by 


r i u=ij)z,  (5.4) 

n . 

v - - *x  . (5.5) 

If  we  set 
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(5.6) 


* - I(5|1)vxu(x)]1/2f(n) 

i 

i 

where  v is  the  kinematic  viscosity,  the  laminar  boundary 

1 i 

layer  equations  reduce  to  the  Falkner-Skan  equation 

2 

f"  + ff"  + B[l-(f')  ] = 0 . (5.7) 

Here  primes  indicate  differentiation  with  respect  to  n 
The  boundary  conditions  for  this  equation  are  as  follows: 

(i)  The  requirement  of  zero  velocity  on  the  plate 
gives: 

^ f (0)  = f ' (0)  = 0 (5.8) 

(ii)  The  requirement  that  the  flow  velocity  approach  the 
free  stream  value  as  z -*•  » implies 

f ' (n)  1 (n  -*■  +00)  • (5.9) 


Properties  of  the  solutions  to  this  boundary  value  problem  have 
been  extensively  investigated  [9]. 

We  solved  (5.7)  - (5.9)  numerically  using  a fourth-order 
Runge-Kutta  shooting  method  together  with  the  mappings  introduced 
earlier.  Some  results  of  our  computations  are  given  in  Table  5* 

In  this  table,  we  give  the  errors  in  the  approximate  solution 
for  f * (n ) at  n = 1 for  the  restricted  domain,  exponential 
mapping  and  algebraic  mapping  for  a variety  of  choices  of  scale 
L for  fixed  total  number  of  grid  points.  It  is  apparent  from  these 
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results  that  the  algebraic  mapping  allows  substantial  economies  in 
the  numerical  solution  of  the  Falkner-Skan  equation.  As  in  Sec.  4, 
the  efficiency  of  the  mapping  methods  is  most  important  in  the 
context  of  multidimensional  numerical  hydrodynamics,  where 
economy  in  resolution  is  vital. 

6.  Wave  Equation 

In  this  section  we  study  the  application  of  mapping  to 
the  problem 

utt  * uxx  (0  < x < »,  t > 0) 

u(0, t)  = f(t)  (t  > 0)  (6.1) 

u(x,0)  = ufc(x,0)  =0  (0  < x < ®)  . 

The  exact  solution  to  this  problem  is 

( f(t-x)  X < t 

u (x , t ) = { (6.2) 

( 0 x > t 

which  represents  an  outgoing  wave  at  x = + « . 

If  the  semi-infinite  interval  0 < x < ® is  replaced 
by  the  finite  interval  0 £ x £ L , boundary  conditions  must 
be  applied  at  x = L . If  the  boundary  condition  that  is  applied 
is  simply 

u (L,  t)  » 0 

waves  reflect  from  the  boundary;  the  exact  soltuion  to  (6.1) 
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on  the  interval  0 <_  x < L with  u(L,t)  =0  is 

00 

u(x,t)  » l [f (t-x-2nL)  - f (t+x-2(n+l)L) ] (6*3) 

n*0 

where  f(s)  = f(s)  if  s > 0 and  f(s)  = 0 if  s £ 0 . The 
solution  (6.3)  consists  of  the  original  outgoing  wave  (6.2)  to- 
gether with  reflected  waves  that  begin  to  appear  at  t = nL, 
n = 1,2,...  # If  x is  near  0 the  solution  (6.3)  is  a good 
approximation  to  the  exact  solution  (6.2)  only  for  t < 2L  when 
the  first  reflected  wave  arrives  at  x = 0 . 

There  are  other  boundary  conditions  that  can  be  applied 
at  x = L that  do  not  yield  reflected  waves.  If  we  set 

ut  + ux  = 0 (6.4) 


■ i 


t 

l ! 

r ! 


at  x = L , the  exact  solution  (6.2)  is  recovered.  However, 
we  do  not  enter  into  a discussion  of  radiation  boundary  con- 
ditions like  (6.4)  here;  a full  discussion  of  them  will  be  given 
elsewhere. 

The  mappings  (2.7),  (2.8)  can  also  be  applied  to  the  wave  prob- 
lem (6.1).  For  example  with  the  algebraic  map  (2.8)  (6.1)  becomes 


32u  _ (1-z)2  3 ,,  _,2  3u 

^ 15  (1-*>  37 


(0  < z < 1,  t > 0) 


(6.5) 


U 


To  analyze  these  results,  we  note  that  there  are  basically 
two  kinds  of  errors  in  the  numerical  solution:  first,  there  are 

local  truncation  errors  in  the  representation  of  the  wave  prop- 
agating through  x ; and,  second,  there  are  the  reflected  waves 
from  the  artificial  boundary  at  z = 1 . In  practice,  the  first 
kind  of  error  is  kept  to  within  a few  percent  with  a second 
order  difference  scheme  so  long  as  the  grid  resolution  Ax 
satisfies 

Ax  < £ (6.6) 

where  X is  the  wavelength  and  M is  a number  of  order  10; 
roughly  speaking,  at  least  10  grid  points  per  wavelength  are 
required  for  an  accurate  calculation. 

The  error  due  to  the  reflected  waves  is  more  serious; 
when  the  first  reflected  wave  arrives  the  solution  is  completely 
wrong!  It  turns  out  in  practice  that  waves  are  reflected 
appreciably  either  from  the  boundary  at  x = L if  there  is 
no  map  or  with  a map,  from  the  points  x at  which 

Ax  * ^X  (6.7) 

When  the  local  grid  spacing  gets  larger  than  roughly  ^ X 
as  it  must  for  z close  to  1 , waves  can  no  longer  be  re- 
solved on  the  mapped  grid.  Since 


X when 


Thus 


with  the  algebraic  map  (assuming  2L  <<  NX)  and  near 


x = L In 


NX 

2L 


(6.9) 


with  the  exponential  map.  These  results  are  consistent  with  the 
numerical  results  plotted  in  Figs.  5-7. 

Now  we  pose  the  following  question.  Suppose  that  N , 
the  number  of  grid  points,  is  fixed  and  that  we  wish  to  calculate 
k wavelengths  of  the  solution  (6.2)  accurately  for  as  long  a 
time  as  possible.  Since  the  scale  of  the  maps  is  a free  para- 
meter, we  can  choose  L to  maximize  the  time  interval  of  accuracy. 
The  question  is:  for  how  long  do  the  various  mappings  maintain 

an  accurate  solution?  We  analyze  three  cases: 

(i)  z = x/L  . In  this  case,  the  calculation  is  accurate 
until  roughly  t = 2L  when  the  first  reflected  wave  arrives 
(assuming  L >>  kX) . Since  we  must  choose  | = Ax  $ X/M 
(where  H » 10)  it  follows  that 


N 
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accuracy 


_ 2XN 
M 


r 


k . 

i 


(ii)  z = x/ (x+L)  . In  this  case  (6.8)  shows  that  the  first 
reflected  wave  arrives  at  x = 0 at  t = ^2NXL  . However, 
small  local  truncation  error  for  the  first  k wavelengths 
about  x = 0 requires  that,  for  z near  0 , Ax  < X/M  . 
Thus,  we  must  also  require  L $ NX/M  . Thus,  the  optimal 
time  of  accuracy  is  roughly 


accuracy 


■VI 


NX 


(iii)  z = 1 


_ a-x/L 


Eq.  (6.9)  implies  that  the  first  reflected 


waves  reach  x = 0 at  t = 2L  Jln(NX/2L).  As  in  (i)  and  (ii)  above, 
L < NX/M  so  that 


t = £ NX  Jin  (M/2)  . 

accuracy  M 


These  results  show  that,  if  high  accuracy  is  desired  so  M >>  1 , 
then  the  algebraic  map  (ii)  gives  the  longest  time  of  accuracy 
followed  by  the  exponential  map  (iii)  and  the  restricted  domain 
(i) . However,  for  the  cases  plotted  in  Figs.  5-7  the  exponential 
map  gives  the  largest  value  of  taccuracy  * 
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7.  Burger's  equation 

In  this  section  we  consider  the  utility  of  mappings  for 
the  problem 


3u  3u  - 32u 


(0  < x < 


t > 0) 


(7.1) 


u(x,0)  = f(x) 
u (0 , t)  = g(t)  . 


(7.2) 

(7.3) 


Because  of  the  nonlinearity  of  this  problem,  there  are  choices 
of  f(x)  and  g(t)  for  which  mapping  is  an  effective  technique 
to  solve  the  problem  numerically  and  other  choices  of  f (x) 
and  g(t)  for  which  mapping  is  not  effective.  In  general, 
if  the  solution  approaches  a constant  as  x -*•  « uniformly  in 
t , mapping  is  effective;  if  the  solution  violates  this  condition 
mapping  is  not  effective.  We  will  now  give  examples  to  illustrate 
these  remarks. 

If  f (x)  = e-x  - 1 , g(t)  = 0 , then  the  solution  to 
(7.1-3)  satisifes 


u(x,t)  <v-  - tanh 


( t -►  -H») 


(7.4) 


For  this  problem,  mapping  methods  work  well.  In  Fig.  8 we  plot 
the  results  of  numerical  integration  of  (7.1-3)  using  a second- 
order  centered  difference  scheme  together  with  the  algebraic 
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map  (2.8)  with  L = 1 . The  results  plotted  in  Fig.  8 illustrate 
the  way  in  which  the  asymptotic  solution  (7.4)  is  achieved  in 
time. 

On  the  other  hand,  if  f(x)  = 0 , g(t)  = 1 , then  the 
solution  to  (7.1-3)  is  given  asymptotically  by  a propagating 
wave 

x-1t 

u(x,t)  - 1 - tanh  (t  -*■  «)  . (7.5) 


For  any  of  the  numerical  methods  with  a fixed  grid  to 
treat  the  boundary  at  x = ® , there  will  be  a time  T beyond 

which  the  spatial  resolution  is  not  sufficient  to  reproduce 
u(x,t)  accurately.  If  the  domain  is  restricted  to  0 <_  x <_  L , 
the  asymptotic  result  (7.5)  implies  that  the  solution  is  accurate 
only  for 

t < 2L  . 

If  one  of  the  mappings  is  used  the  solution  becomes  in- 
accurate when  the  shock  location  at  x = ^-t  is  within  a 
region  of  poor  x resolution.  On  the  basis  of  the  results 
in  Sec.  6 [see  (6.7)]  we  expect  that  large  errors  will  result 
when  the  shock  thickness  2v  is  larger  than  Ax  . For 
both  the  algebraic  and  exponential  maps 


- 


V- 


JEl 


dx 

3i 


(1-z) 


where  a = 1 for  the  exponential  map  and  a = 2 for  the 
algebraic  map.  Therefore  if  Az  = 1/N  , the  effective 
resolution  Ax  = 2v  when 


1 - 2 * (2Sv> 


1/a 


where  we  assume  that  L <<  2Nv.  Thus,  the  solution  using 
the  algebraic  map  is  expected  to  deteriorate  at 


alg 


2Nv 


(7.6) 


I. 

* 

r>  i 

r 


while  the  exponential  map  should  give  results  that  deteriorate 
near 


xexp  * 2Nv  (2^>  TT 


(7.7) 


* 


« i 

I 

- 

% 


Note  that  (7.6-7)  shows  the  explicit  dependence  of  xa^g  and 

x^ on  the  small  parameter  L/2Nv  ; observe  that  x <<  x . 

exp  exp  alg 

These  predictions  are  consistent  with  the  results  plotted  in 
Figs.  9-10. 


25 


— 


8.  Conclusion 


We  conclude  from  the  examples  given  above  that  mappings 
are  an  effective  way  to  solve  problems  in  infinite  domains 
provided  that  the  solution  is  simple  at  infinity.  If  the 
solution  oscillates  as  x -►  00  then  00  must  be  an  essential 
singularity  of  the  solution  and  mappings  fail.  Unfortunately 
this  implies  that  mappings  are  nearly  useless  for  many  important 
physical  problems. 

When  mapping  is  applicable,  the  proper  choice  of  mapping 
should  be  based  on  the  criterion  that  the  solution  to  the 
problem  be  smooth  in  the  mapped  coordinate.  For  many  problems 
this  criterion  favors  algebraic  mappings  over  exponential  mappings. 


Table  1.  Errors  in  the  numerical  solution  of  the  Rayleigh  shear  flow  at  x = 1.0  as 
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Table  2.  Convergence  of  approximations  to  the  eigenvalue 
X = 4.5  of  the  quantum-mechanical  harmonic 
oscillator  using  N+l  Chebyshev  polynomials. 


Map 

L 

N 

X 

Truncation 

4.0 

10 

5.20628  14857 

4.0 

30 

5.20628  12800 

-L  < x £ L 

8.0 

10 

4.57205  38006 

8.0 

20 

4.50000  00394 

8.0 

30 

4.50000  00395 

16.0 

20 

4.56679  36320 

16.0 

30 

4.50000  17110 

32.0 

30 

5.47165  94003 

Algebraic 

4.0 

20 

4.50029  00880 

2 

4.0 

30 

4.49999  99641 

z - 2 J - 1 

x^  + L 

8.0 

10 

4.56858  45536 

8.0 

20 

4.50002  73879 

8.0 

30 

4.49999  99985 

16.0 

20 

4.50000  02905 

16.0 

30 

4.50000  00000 

32.0 

10 

4.49656  47888 

32.0 

20 

4.49999  98895 

32.0 

30 

4.50000  00000 

64.0 

20 

4.49999  99898 

128.0 

20 

4.49999  96979 

256.0 

20 

4.62092  09932 

I 


} 


Table  3.  Values  of  the  eigenvalues  X of  the  quantum' 
. mechanical  harmonic  oscillator  obtained  with 
the  algebraic  map  and  N = 20  (21  Chebyshev 
polynomials) . 
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.5000  0171 

2.5016  8050 

4.4889  8957 

2 

.5000  0004 

2.5001  2172 

4.4968  0153 

4 

.5000  0000 

2.5000  0347 

4.5002  9009 

8 

.5000  0000 

2.5000  0000 

4.5000  2739 

16 

.5000  0000 

2.5000  0000 

4.5000  0029 

32 

.5000  0000 

2.5000  0000 

4.4999  9989 

64 

.5000  0000 

2.5000  0000 

4.4999  9999 

128 

.5000  0000 

2.5000  0000 

4.4999  9970 

256 
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Table  5.  Errors  in  solution  of  the  Falkner-Skan  equation 
using  11  grid  points. 


Method 

Z. 

Error  in  f ' (1) 

e = o 

£'  (1)  = 0.4606  3259 

Restricted  Domain 

1.0 

5.4  x lo"1 

2.5 

5.2  x lo"2 

5.0  ' 

1.3  x io"4 

7.5 

• 

10. 

* 

Exponential  Map 

1.0 

-1.1  x io-1 

2.5 

-1.2  x io"1 

5.0 

-1.0  x io"1 

7.5 

-8.2  x io"2 

10. 

-6.9  x io’2 

Algebraic  Map 

1.0 

5.8  x 1<T6 

2.5 

-4.7  x io“7 

5.0 

-2.8  x io”6 

7.5 

-4.0  x IO”6 

10. 

-6.7  x io”5 

6 -*  -.1 

f ' Cl)  * 0.3630  0791 

Restricted  Domain 

5.0 

4.2  x io”4 

7.5 

• 

Exponential  Map 

5.0 

-8.8  x io"2 

7.5 

-7.3  x io”2 

Algebraic  Map 

5.0 

-4.6  x IO-6 

7.5 

-6.4  x 10”6 

6 - .1 

f’  (1)  = 0.5274  2343 

Restricted  Domain 

2.5 

3.2  x 10”2 

5.0 

1.4  x 10“4 

7.5 

• 

Exponential  Map 

2.5 

-1.3  x 10”6 

5.0 

-1.1  x io”1 

7.5 

-8.7  x IO-2 

•Is 

Algebraic  Map 

2.5 

-5.5  x 10  ° 

5.0 

-1.1  x 10“5 

7.5 

-1.6  x IO*5 

* No  acceptable  (monotonic)  solution. 


Figure  1.  Coordinate  system  for  flow  past  a semi-infinite  flat  plate 


Figure  2.  Variation  of  z versus  x for  the  exponential  map  with  various 
values  of  L.  The  points  on  the  curves  indicate  values  of  x 


Variation  of  z versus  x for  the  algebraic  map  with  various  values 


Cbebyshev  approximation  (3.3)  with  N = 18  is  not  distinguishable  from  the  exact 
infunction  u(x)  = -(1  - x2)  exp(-x2/4) . 


„ r. 


Figure  5.  Plot  of  the  solution  to  the  wave  equation  (6.1)  at  x = 1 as  a function  of 
t using  the  restricted  domain  method  with  L = 5. 


Figure  6.  Plot  of  the  solution  to  the  wave  equation  (6.1)  at  x = 1.25  as  a function  of 
t using  the  algebraic  map  with  L = 5. 


Plot  of  the  solution  to  the  wave  equation  (6.1)  at  x = (-5  In  8)/1.15  as 
function  of  t using  the  exponential  map  with  L = 5. 
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